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Abstract 

This  paper  addresses  the  problem  of  nonlinear  multivariate  root  finding.  In  an  earlier  paper  we  describe 
a  system  called  Newton  which  finds  roots  of  systems  of  nonlinear  equations  using  refinements  of  interval 
methods.  The  refinements  are  inspired  by  AI  constraint  propagation  techniques.  Newton  is  competitive 
with  continuation  methods  on  most  benchmarks  and  can  handle  a  variety  of  cases  that  are  infeasible  for 
continuation  methods.  This  paper  presents  three  “cuts”  which  we  believe  capture  the  essential  theoretical 
ideas  behind  the  success  of  Newton.  This  paper  describes  the  cuts  in  a  concise  and  abstract  manner  which, 
we  believe,  makes  the  theoretical  content  of  our  work  more  apparent.  Any  implementation  will  need  to 
adopt  some  heuristic  control  mechanism.  Heuristic  control  of  the  cuts  is  only  briefly  discussed  here. 
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1  Introduction 

In  this  paper  we  address  the  [)roblem  of  finding  solutions 
to  large  systems  of  nonlinear  equations.  This  is  an  old 
problem  with  many  applications  and  a  large  literature. 
In  engineering  applications  it  is  generally  sufficient  to 
find  a  floating  point  value  for  each  variable  such  that  the 
given  equations  are  satisfied  up  to  the  accuracy  of  float¬ 
ing  point  computations.  The  problem  of  determining 
whether  such  a  floating  point  solution  exists  is  NP  com¬ 
plete.  Our  work  emphasizes  refinements  of  constraint 
propagation  techniques  for  NP  complete  problems  that 
have  been  developed  in  the  AI  and  logic  programming 
communities  [1,  2,  3]. 

In  an  earlier  paper  we  describe  an  implemented  sys¬ 
tem  called  Newton  [4].  Newton  uses  a  branch  and  prop¬ 
agate  algorithm  whose  propagation  phase  manipulates 
upper  and  lower  bounds  on  the  variables  appearing  in  the 
given  equations.  Algorithms  which  manipulate  upper 
and  lower  bounds  are  called  interval  methods  —  Newton 
is  based  on  interval  methods  [5].  A  well-studied  alter¬ 
native  to  interval  methods  is  the  so-called  “continuation 
method”  or  “homotopy  technique”  [6].  Table  1  summa¬ 
rizes  our  previously  published  results  comparing  Newton 
with  earlier  interval  based  systems  and  with  continuation 
methods.^  The  benchmarks  were  taken  from  papers  on 
numerical  analysis  [12],  interval  analysis  [8,  9,  11],  and 
continuation  methods  [14,  6,  13,  10].  The  table  shows 
that  Newton  outperforms  an  earlier  interval  based  sys¬ 
tem  and  is  competitive  with  continuation  methods  on 
the  benchmarks  tried.  The  novelty  of  Newton  resides 
in  the  use  of  bound  propagation  techniques  which  are 
stronger  than  those  used  in  previous  interval  systems. 
Newton  can  solve  a  variety  of  problems  for  which  con¬ 
tinuation  methods  are  infeasible  (those  benchmarks  with 
total  degree  greater  than,  say,  2^^). 

Our  earlier  paper  on  Newton  describes  the  algorithm 
used  in  sufficient  detail  to  allow  replication  of  the  empiri¬ 
cal  results  [4].  This  includes  a  detailed  description  of  the 
heuristics  used.  The  details  tend  to  obscure  what  we  feel 
are  the  essential  technical  ideas  behind  the  algorithms. 
In  this  paper  we  give  an  abstract  presentation  of  three 

^The  Newton  system  was  run  on  a  Sun  Sparc  10  work¬ 
station  to  obtain  all  solutions  to  each  benchmark  system  of 
equations.  The  column  labeled  Newton  gives  the  run  time  (in 
seconds)  of  the  newton  system  on  a  given  benchmark  prob¬ 
lem.  The  column  labeled  HRB  gives  timings  for  a  traditional 
interval  method  using  Hansen-Segupta’s  operator,  range  test¬ 
ing,  and  branching.  This  method  uses  the  same  implemen¬ 
tation  technology  as  ours,  C  code  running  on  a  SPARC  10. 
Some  interval  methods  such  as  [7]  have  better  convergence 
properties  near  solutions.  Our  main  contribution  aims  at 
speeding  up  the  computation  when  far  from  a  solution  and 
hence  comparing  it  to  HRB  is  meaningful.  The  column  labeled 
CONT  gives  timings  for  a  state  of  the  art  continuation  method 
[14]  running  on  a  DEC  5000/200  which  is  perhaps  slightly 
slower  than  a  SPARC  10.  For  each  benchmark,  we  give  the 
number  of  variables  (n),  the  total  degree  of  the  system  (d), 
the  initial  range  for  the  variables,  and  the  run  time  of  each 
system  in  seconds.  A  space  in  a  column  means  that  the  re¬ 
sult  is  not  available  for  the  method.  A  question  mark  means 
that  the  method  does  not  terminate  in  a  reasonable  time  (> 

1  hour).  Further  details  on  these  timings  can  be  found  in  [4]. 
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Table  1:  Summary  of  the  Experimental  Results 


‘"cuts”  where  each  cut  specifies  a  way  of  inferring  a  new 
upper  or  lower  bound  on  a  variable.  These  abstract  cuts, 
while  insufficient  in  themselves  to  allow  exact  replication 
of  our  results,  seem  to  capture  the  essence  of  Newton’s 
ability  to  prune  the  search  space  and  a  description  of 
the  cuts  should  allow  others  to  replicate  the  qualitative 
performance  of  the  Newton  system. 

The  heart  of  Newton’s  algorithm  is  a  propagation  pro¬ 
cess  which  iteratively  improves  known  bounds  on  vari¬ 
ables  using  a  set  of  inference  rules  we  call  interval  cuts. 
An  interval  cut  is  a  method  of  inferring  an  inequality 
of  the  form  x  <  b  ov  x  >  b  where  6  is  a  floating  point 
number.  The  next  section  specifies  three  abstract  cuts 
used  in  the  Newton  system.  Section  3  briefly  outlines 
how  these  cuts  can  be  used  in  a  branch  and  propagate 
procedure  for  findings  all  roots  of  a  system  of  nonlinear 
equations.  Section  3  is  a  simplification  of  the  algorithm 
actually  used  in  Newton. 

2  Three  Cuts 

In  order  to  define  our  cuts  some  preliminary  terminol¬ 
ogy  is  needed.  Our  algorithms  work  with  floating  point 
numbers  wffiich  we  assume  to  be  a  finite  set  of  ratio¬ 
nal  numbers  extended  with  both  a  positive  and  nega¬ 
tive  infinity,  denoted  Too  and  — cxo  respectively.  We 
will  generally  use  the  term  “interval”  to  refer  a  closed 
interval  of  the  form  [a,  6]  where  a  and  b  are  floating 
point  numbers.  Our  algorithms  also  work  with  symbolic 
expressions  constructed  from  floating  point  constants, 
symbolic  variables,  addition,  subtraction,  multiplication, 
and  “power  expressions”  of  the  form  where  e  is  an 
expression  and  n  is  a  positive  integer  constant.  Power 


expressions  improve  the  accuracy  of  interval  evaluation. 
We  will  generally  use  the  term  “expression”  to  refer  to 
symbolic  expressions  of  this  form.  We  will  use  the  term 
“box”  to  refer  to  an  assignment  of  intervals  to  symbolic 
variables.  For  box  B  and  variable  x  we  write  B{x)  to  de¬ 
note  the  interval  assigned  to  x  in  the  box  B,  A  “point** 
is  an  assignment  of  particular  floating  point  numbers  to 
each  variable.  Clearly,  if  we  are  working  with  n  variables 
then  a  box  is  a  rectilinear  subset  of  .  We  allow  B{x) 
to  be  interval  of  the  form  [a,  a].  A  point  can  be  viewed 
as  a  box  in  which  all  intervals  are  of  this  form. 

All  interval  methods  are  based  on  interval  evalua¬ 
tion  —  an  evaluation  process  which  computes  an  inter¬ 
val  value  for  an  expression  from  given  intervals  for  the 
variables  it  contains.  The  most  straightforward  interval 
evaluation  method  is  obtained  by  associating  each  of  the 
operations  of  addition,  subtraction,  multiplication,  and 
powers  with  a  corresponding  operation  from  intervals  to 
intervals.  More  specifically  we  have  the  following  opera¬ 
tions  on  intervals. 


j  NE(e.  B{x:=(c  ,c’]]) 


Line  of  slope  Du  Line  of  Slope  Dl 


Figure  1:  A  Newton  cut  deriving  x  <  c'  — 

procedure  can  also  be  viewed  as  an  inference  rule  for  de¬ 
riving  new  bounds  independent  of  heuristic  methods  for 
finding  useful  inferences. 


[a,  b]  -(-  [c,  d]-[a-\-c,  b d] 

[a,  6]  —  [c,  d]  ==  [a  —  d,  b  —  c] 

[a,  b]  ♦  [c,  d]  =  [min{ac,  ad,  be,  bd),  max{ac,  ad,  be,  bd)] 


[a,  br  = 


[a^,  6”]  n  odd 

’0,mai(a”,  6”)]  n  even;  0  6  [a,  6] 

’min{a^ ,b^),  max(a”,6”)]  n  even;  0  ^  [a,  b] 


We  use  the  notation  NE{e,  B)  to  denote  the  inter¬ 
val  value  of  an  expression  e  when  variables  are  assigned 
the  interval  values  given  in  B  and  operations  are  inter¬ 
preted  as  operations  on  intervals  according  to  the  above 
rules  (lYE  stands  for  “natural  evaluation”).  For  exam¬ 
ple  if  B{x)  is  [—1,  1]  then  N E{x  +  x,  B)  is  the  inter¬ 
val  [-2,  2];  NE{x*x,  B)  is  the  interval  [—1,  1];  but 
NE{x^,  B)  is  [0,  1].  For  most  nontrivial  expressions 
the  interval  NE{e,  B)  is  larger  than  the  actual  range  of 
values  of  e  over  points  in  B.  If  B{x)  is  [—1,  1]  we  have 
that  X E{i'“  —  X,  5)  is  [—1.  2]  while  the  actual  range  of 
values  in  this  box  is  [—1/4,  2].  In  our  implementation 
outward  rounding  is  used  in  the  above  rules  —  in  arith¬ 
metic  for  computing  lower  bounds  we  round  down  and 
in  arithmetic  for  computing  upper  bounds  we  round  up. 
This  provides  a  guarantee  that  NE{€,  B)  contains  the 
actual  range  of  values  of  e  independent  of  numerical  er¬ 
rors.  This  guarantee  is  often  expressed  as  the  statement 
that  interval  evaluation  is  “conservative” . 

We  will  assume  a  given  set  of  constraints  of  the  form 
e  <  0  where  e  is  an  expression  as  defined  above.  An 
equation  ci  =  62  can  obviously  be  represented  with  two 
constraints  of  the  form  ei  —  62  <  0  and  62  —  ei  <  0.  For 
technical  reasons  it  will  be  easier  to  work  with  inequali¬ 
ties. 

An  interval  cut  is  a  procedure  which  operates  on  a 
set  of  constraints  and  a  current  box.  An  interval  cut 
reduces  the  given  box  by  deriving  a  new  bound  on  one 
of  the  variables.  Cuts  are  presented  below  as  nondeter- 
ministic  procedures  —  performing  a  cut  often  involves 
guessing  one  or  more  values.  In  practice  these  guesses 
are  generated  by  deterministic  heuristic  methods  such 
as  those  described  in  section  3.  A  nondeterministic  cut 


2.1  Newton  Cuts 

A  Newton  cut  is  represented  schematically  in  figure  1. 
The  notation  used  in  the  figure  is  defined  rigorously  in 
the  following  description.  In  the  Newton  cut,  and  in 
later  cuts,  we  abuse  notation  and  use  de/dx  to  denote 
the  expression  that  is  the  symbolic  derivative  of  e  with 
respect  to  x. 

Newton  Cut:  Let  e  <  0  be  a  constraint,  B  a 
given  box,  x  a  variable  appearing  in  e,  and  let 
[a,  i]  be  the  interval  B{x).  Let  c  and  c'  be 
two  cut  values  such  that  c  <b  and  c'  in  [c,  b]. 

Let  s  be  the  lower  bound  of  NE{e,  B[x  := 
[c',c']]).  We  require  5  >  0.  Let  [Di,  Du]  be 
the  interval  NE{de/dx,  B[x  :=  [c,  6]]).  To 
ensure  there  is  no  solution  with  x  6  [c',6]  we 
require  that  either  Di  >  0  oi  Di  <  0  and 
c'  —  >  6.  If  Dti  <  0  we  can  infer  x  <  c  and 

if  Du  >  0  we  can  infer  x  <  max(c,  d  —  ;^)- 

The  lines  of  slope  Du  and  Di  in  the  figure  1  represent 
the  fastest  possible  rate  of  descent  of  the  value  of  the 
expression  e  as  a  function  of  x  based  on  the  interval  eval¬ 
uation  of  the  derivative  dejdx.  It  should  be  clear  that 
all  values  of  x  inside  the  triangle  defined  by  these  lines 
must  violate  the  constraint  e  <  0  and  can  be  pruned. 

An  interesting  special  case  is  when  c  =  a  and  c'  =  6. 
Whenever  the  lower  bound  of  NE{e,  B[x  :=  6,6]])  is 
greater  than  zero  the  upper  bound  6  can  be  ruled  out  as 
a  value  of  x.  If  we  are  working  with  a  finite  box,  i.e., 
one  where  all  variable  intervals  are  finite,  and  ignoring 
roundoff  errors,  Du  must  be  finite  and  the  choice  of  c  =  a 
and  d  —  b  always  achieves  some  reduction  of  the  upper 
bound.  However,  larger  cuts  are  usually  achieved  by 
other  choices  of  c  and  c'.  Increasing  the  value  of  c  reduces 
the  size  of  the  interval  [c,  6]  and  hence  reduces  the  range 
of  possible  derivative  values  NE{deldx,  B[x  [c,  6]]). 
This  often  gives  a  smaller  value  for  Du  and  hence  a  larger 
cut. 

Decreasing  the  value  of  c'  away  from  6  can  also  give 
a  larger  cut.  As  in  the  statement  of  the  cut,  let  s(c')  be 


the  lower  bound  of  NE{e,  B[x  :=  ir\  c']]).  as  a  function 
of  In  practice  a  reduction  in  the  upper  bound  can 
only  be  achieved  when  s{b)  >  0.  i.e..  the  value  b  ciin  be 
ruled  out.  As  c'  is  reduced  from  6,  the  lower  bound  s(c') 
will  also  typically  be  reduced.  However,  it  typically  falls 
more  slowly  than  Du,  which  is  the  largest  possible  value 
of  de/dx  in  the  interval  [c,  6].  If  ds/dc'  is  smaller  than 
Du  then  reducing  c'  away  from  b  will  result  in  a  larger 
cut.  However  if  we  reduce  E  too  aggressively  then  the  cut 
may  fail  altogether  either  because  we  fail  to  get  s{c')  >  0 
or  we  fail  to  rule  out  feasible  points  with  x  G  b]. 

The  above  cut  allows  for  the  reduction  of  upper 
bounds.  We  call  it  an  upper  bound  cut.  Every  upper 
bound  cut  has  a  dual  lower  bound  cut.  We  will  only 
present  upper  bound  cuts. 

2.2  Snake  Cuts 

Consider  a  given  constraint  e  <  0.  a  given  variable  x 
appearing  in  the  constraint,  a  given  box  5,  and  let  [a,  b] 
be  the  interval  B{x),  Now  consider  a  value  c'  in  the 
interval  [a,  6]  and  consider  the  lower  bound  s  of  the 
interval  NE{e,  B[x  [c',c']])  as  a  function  of  c'.  We 
will  write  s(c')  to  express  s  as  a  function  of  c'  for  fixed 
values  of  e  and  B.  The  snake  cut  is  identical  to  the 
Newton  cut  except  that  de/dx  is  replaced  by  ds/dE? 
It  turns  out  that  the  range  of  values  of  dsjdd  over  the 
interval  [c,  b]  is  often  significantly  less  than  the  range  of 
possible  values  of  values  of  de/dx.  This  fact  allows  cuts 
to  be  larger. 

In  order  to  replace  de/dx  by  ds/dc'  we  need  to  first 
represent  the  function  s(c')  as  a  symbolic  expression  s{x) 
involving  only  the  variable  x.  To  do  this  we  first  sym¬ 
bolically  rewrite  e  as  a  polynomial  P{x)  of  the  form 
.So  +  six  -h  S2X^  +  SnX^  wherc  each  coefficient  Si 
is  an  expression  not  involving  x.  We  then  define  two 
polynomials  in  a?,  P^{x)  and  each  of  which  has 

constant  coefficients  (and  hence  no  variables  other  than 
x).  P~{x)  gives  a  lower  bound  on  e  for  negative  values  of 
X  and  P'^{x)  gives  a  lower  bound  on  e  for  positive  values 
of  X.  For  odd  powers  of  x  the  coefficients  of  P~{x)  are 
taken  to  be  the  upper  bound  of  NE(si,  B)  where  Si  is 
the  corresponding  symbolic  coefficient  of  the  polynomial 
P{x).  In  all  other  cases  the  coefficients  of  P~{x)  and 
P'^{x)  are  taken  to  be  the  lower  bound  of  NE(si^  B). 
A  function  5(a:)  can  now  be  represented  symbolically  by 
the  conditional  expression  if  (x  >  0.  P+(x),  P-(x)) 
which  has  the  property  that  it  equals  P'^{x)  for  posi¬ 
tive  values  of  x  and  P~{x)  for  other  values  of  x.  The 
interval  evaluation  function  NE  and  the  symbolic  differ¬ 
entiation  operators  can  be  easily  extended  to  handle  con¬ 
ditional  expressions.  The  symbolic  expression  s{x)  does 
not  correspond  exactly  to  the  function  s(c')  as  defined 
above  because  rewriting  an  expression  as  a  polynomial 
in  X  changes  the  interval  returned  by  interval  evaluation. 
However,  the  symbolic  expression  s(x)  does  give  a  valid 
lower  bound  on  the  values  of  e  for  points  in  the  box  B 
as  a  function  of  the  value  for  x. 


^In  informal  discussions  we  began  calling  the  function 
.s(c')  a  “snake”  and  the  cut  based  on  an  analysis  of  this  func¬ 
tion  became  known  as  the  snake  cut. 


Snake  Cut:  Let  e  <  0,  x,  P,  [a,  b]  and  the 
symbolic  expression  s(x)  be  defined  as  above. 

Let  c  and  P  be  two  additional  parameters  as 
in  the  Newton  cut  with  c  <  b  and  P  in  [c,  6]. 

We  require  that  s(c')  >  0.  Let  [Di,  Du]  be 
the  interval  NE{ds/dx,  B[x  :=  [c,  6]]).  To 
ensure  there  is  no  solution  with  x  E  [P ,  b]  we 
require  that  either  Di  >  0  or  Di  <  0  and 
P  -  >  b.  If  Du  <  0  we  can  infer  x  <  c  and 

if  Du  >  0  we  can  infer  x  <  max(c,  P  —  ;^). 

To  see  the  advantage  of  the  snake  cut  over  the  Newton 
cut  consider  the  expression  zx  —  2.  <  0  where  the  box  B 
assigns  both  x  and  r  the  interval  [1,  4].  It  is  not  difficult 
to  see  that  since  z  is  in  [1,  4]  and  zx  <  2  we  have  x  <  2. 
We  compare  a  Newton  cut  with  a  snake  cut  where  we 
take  c  ==  1  and  c'  ~  4  in  both  cuts  (i.e,,  c  =  a  and 
P  =  b).  In  the  Newton  cut  we  have  that  NE{de/dx,  B) 
equals  NE[z,  B)  which  is  [1,  4].  So  the  infered  bound 
isx<4-|  =  3|.  However  the  lower  bound  polynomial 
s{x)  (for  positive  values  of  x)  is  x  — 2.  So  NE{ds/dx,  B) 
equals  NE{1,  B)  equals  [1, 1]  and  the  infered  bound  is 
X  <  4  —  I  =  2,  an  optimal  cut. 

Unfortunately  snake  cuts  require  rewriting  the  con¬ 
straint  in  a  way  that  makes  interval  evaluation  less  ac¬ 
curate.  In  our  experience  snake  cuts  are  most  useful  in 
very  large  boxes  and  Newton  cuts  become  more  useful 
as  the  boxes  get  smaller. 

2.3  Combination  Cuts 

The  combination  cut  is  used  to  gain  some  of  the  power  of 
multivariate  Newton's  method  once  the  search  has  begun 
to  converge  on  a  solution.  Interval  versions  of  multivari¬ 
ate  Newton’s  method  are  widely  used  in  interval  systems, 
e-g-,  [8,  7].  The  Newton  cuts  we  use  involve  two  ideas. 
First,  we  compute  combined  constraints  by  inverting  the 
“center  Jacobian”  matrix  of  the  constraint  system.  Sec¬ 
ond,  in  the  combination  cuts  we  use  a  different  method 
of  interval  evaluation  called  Taylor  evaluation  which  is 
more  accurate  for  small  boxes. 

Like  the  natural  evaluation  operator  NE,  Taylor  in¬ 
terval  evaluation  takes  an  expression  e  and  a  box  B  and 
returns  an  interval  containing  the  range  of  values  of  e 
on  points  in  B.  Taylor  interval  evaluation  tends  to  be 
more  accurate  in  the  case  where  B  is  small.  To  define 
the  Taylor  interval  evaluation  we  first  define  Be  to  be  the 
center  of  the  box  B.  i.e.,  Bc{x)  is  the  interval  [c,  c]  where 
c  is  the  center  of  the  interval  B{x).  We  use  the  notation 
Xi  e  e  to  indicate  that  Xi  is  a  variable  appearing  in  the 
expression  e. 

T£(€,  B)  =  NE{e,  Be)  +  ^  NE^dejdxi^B)  *{,B{xi)  - 

All  arithmetic  operations  in  the  above  expression  are 
operations  on  intervals.  Taylor  interval  evaluation  is 
more  accurate  than  natural  interval  evaluation  in  small 
intervals.  For  example,  consider  the  polynomial  x^  —  x 
where  B{x)  is  the  interval  [xc  —  e,  x^  +  e]  where  Xc>  \ 
and  e  is  a  small  positive  number.  We  have 

NE{x\-xc.  J3)  =  [(xc -Jrc)~(2a;c  +  l)e+e^,  (x^  -  xc)  +  (2xc  +  l)e+ 
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TB(xl~Xc,  J5)  =  {(j'-j:c')-.:2xc-1)€+2€^.  (x^  -  Xc  )  +  (2a;c  -  1 )«  +  2€“] 

For  small  values  of  e  the  width  of  the  first  interval  is 
approximately  2(2xc  +  l)e  while  the  width  of  the  second 
interval  is  approximately  2{2xc  —  l)e.  Since  we  have  as¬ 
sumed  Xc  >  ^  the  first  interval  is  wider  than  the  second. 
(If  we  assume  Xc  <  k  then  the  interval  product  in  the 
Taylor  interval  calculation  gives  a  different  symbolic  an¬ 
swer  and  the  Taylor  interval  is  still  smaller.)  The  natu¬ 
ral  evaluation  does  not  take  into  account  the  correlations 
between  the  two  terms  in  the  polynomial  x“  —  x.  More 
generally,  if  P{x)  is  a  polynomial  in  x  and  B{x)  is  an  in¬ 
terval  of  the  form  [x^-e,  Xc4-e]  then  to  first  order  in  e  we 
have  that  the  width  oiTE(P[x),  e)  is  2P\xc)€  where  P' 
is  the  symbolic  derivative  of  P.  P\xc)  is  nonzero  then  as 
e  goes  to  0  the  ratio  of  the  width  of  TE{P{x),  B)  to  the 
width  of  the  true  range  of  values  of  P{x)  approaches  1 
—  Taylor  interval  evaluation  becomes  exactly  accurate 
in  the  limit  of  small  intervals.  As  the  above  example 
shows,  this  is  not  true  of  the  natural  evaluation. 

Combination  cuts  build  new  constraints  which  are  lin¬ 
ear  combinations  of  given  constraints.  If  we  are  given 
ei  <  0  and  <  0  then  for  any  positive  a  and  we  can 
infer  aei  jSe^  <  0.  A  combination  cut  is  identical  to 
a  Newton  cut  except  that  it  uses  a  constraint  derived  as 
a  linear  combination  of  input  constraints  and  uses  TE 
rather  than  NE  in  the  computation  of  s. 

Combination  Cut:  Let  x  be  a  variable,  B  a 
box,  and  [a,  6]  the  interval  B{x),  Let  c  and  c' 
be  two  additional  parameters  with  c  <  b  and 
c'  in  [c,  6].  Let  ei  <  0,  £2^0?  •  •  •  i  en  <  0  be 
a  set  of  constraints  such  that  each  contains 
X.  Let  u  be  a  linear  combination  Q:iei-l-02e24- 
. .  where  each  a*  is  positive.  Let  s  be 

the  lower  bound  of  TE{u^  B[x  :=  [P,  c']]). 

We  require  s  >  0.  Let  [D/,  Du]  be  the  in¬ 
terval  NE{du/dx,  B).  To  ensure  there  is  no 
solution  with  x  G  [c^6]  we  require  that  either 
Di  >  0  or  Di  <  0  and  c'  —  >  b.  If  Du  <  0 

we  can  infer  x  <  c  and  if  Du  >  0  we  can  infer 
X  <  max{c,  d  —  ^). 

To  see  the  power  of  the  combination  cut  consider  the 
two  constraints  x  -f  y  —  2  <  0  and  x  —  y  <  0.  If  both 
B{x)  and  B{y)  are  the  interval  [0,  2]  then  no  Newton 
or  snake  cut  can  be  used  to  improve  the  bounds  —  each 
constraint  is  “satisfied  at  the  end  points” .  For  example 
NE{x  +  2/  -  2,  B[x  :=  [2,  2]])  is  [0,  2]  so  x  =  2  can  not  be 
ruled  out  on  the  basis  of  this  constraint  alone.  To  reduce 
the  box  B  we  must  examine  both  constraints  simulta¬ 
neously.  Adding  the  two  constraints  together  (with  unit 
coefficients)  gives  2x^2  <  0.  This  immediately  gives  the 
tighter  bound  x  <  1.  If  we  are  given  the  four  inequality 
constraints  representing  x  -f  y  =  2  and  x  —  y  —  0  then 
four  combination  cuts  can  be  used  to  solve  for  x  and  y 
exactly.  A  method  for  finding  appropriate  combination 
cuts  is  described  in  the  next  section. 

3  Branch  and  Propagate  Root  Finding 

To  find  a  feasible  point  for  a  set  of  inequality  constraints, 
or  a  root  to  a  set  of  equations,  one  can  use  a  simple 


branch  and  propagate  procedure.  One  starts  with  a  box 
large  enough  to  contain  all  points  of  interest.  One  then 
reduces  the  box  by  applying  cuts  to  derive  new  bounds. 
At  some  point  one  may  decide  to  branch.  To  branch  one 
selects  a  variable  x.  Let  [a, ’6]  be  the  interval  B{x)  where 
B  is  the  box  in  use  at  the  time  of  the  branch.  One  then 
“splits”  the  box  into  the  two  branches  B[x  :=  [a, 
and  B[x  [^^,6]].  One  then  recursively  computes  a 
solution  inside  the  first  box,  or  if  there  is  no  such  solu¬ 
tion,  a  solution  inside  the  second  box.  The  procedure 
terminates  when  the  current  box  is  sufficiently  small, 
e.g.,  every  variable  is  assigned  an  interval  of  width  less 
than,  say,  10“^.  Newton  uses  a  simple  round-robin  strat¬ 
egy  to  select  the  variable  to  split  —  down  any  branch  of 
the  search  tree  Newton  splits  a  variable  that  has  gone 
the  longest  time  without  splitting. 

The  basic  branch  and  propagate  procedure  can  be 
used  with  a  wide  variety  of  heuristics  for  finding  useful 
cuts  and  for  determining  when  to  branch.  Our  experi¬ 
ence  with  Newton  indicates  that  one  useful  heuristic  is 
to  only  perform  cuts  that  result  in  at  least  a  10%  re¬ 
duction  of  an  interval.  This  ensures  that  the  number  of 
cuts  in  a  single  propagation  phase  is  at  worst  linear  in 
the  number  of  variables  in  the  constraints.  However,  this 
same  heuristic  will  force  the  procedure  to  branch  before 
all  possible  cuts  have  been  exhausted.  In  our  experience 
such  “early  branching”  is  usually  desirable. 

Now  consider  selecting  the  values  of  c  and  c'  in  a  cut 
on  variable  x  with  interval  [a,  6].  A  simple  strategy  is 
to  select  c  to  be  6  —  for  some  non-negative  integer  n 
and  to  select  d  to  be  the  midpoint  of  [c,6].  In  practice 
we  need  only  consider  values  of  n  up  to  3  since  we  are 
only  interested  in  cuts  which  reduce  the  interval  by  at 
least  10%,  For  a  given  variable  we  can  start  with  n  =  0 
and  if  that  cut  fails  increase  n  (up  to  3)  looking  for  a 
viable  cut. 

Combination  cuts  require  selecting  not  only  c  and  d 
but  also  a  linear  combination  of  the  constraints.  Be¬ 
cause  computing  appropriate  coefficients  for  the  linear 
combination  can  be  expensive  we  suggest  performing  all 
possible  Newton  and  snake  cuts  before  attempting  com¬ 
bination  cuts.  In  performing  combination  cuts  we  as¬ 
sume  that  the  input  is  given  as  a  set  of  equations  of  the 
form  €(  =  0.  This  allows  one  to  compute  combined  con¬ 
straints  of  the  form  =  0  where  each  is  a  linear 
combination  of  the  and  where  is  intended  to  con¬ 
strain  X  independent  of  the  value  of  other  variables.  To 
compute  the  combined  constraints  u^;  let  B  be  the  box 
existing  at  the  time  we  are  computing  the  combined  con¬ 
straints.  We  define  the  “center  Jacobian”  matrix  [Dtj] 
by 

Bij  =  the  midpoint  of  the  interval  NE{dei/dxj,  B). 

Assuming  that  the  matrix  [Bij]  is  nonsingular  one  can 
compute  [Aij]  as  of  [Bij].  We  then  define  the  expression 
«r,  by 

n 

k  =  \ 

Inside  the  box  D,  and  especially  for  reasonably  small 
boxes,  we  have  that  duxjdxj  is  near  1  if  z  =  j  and  near 
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0  otherwise.  So  Ux  =  0  is  primarilly  a  constraint  on 
X.  Once  the  combined  constraints  of  the  form  Ux  =  0 
have  been  computed  we  can  perform  combination  cuts 
on  this  fixed  set  of  combined  constraints  selecting  c  and 
as  specified  above.  In  our  experience  it  is  best  to 
perform  only  a  single  matrix  inversion  in  any  one  propa¬ 
gation  phase  of  the  procedure  —  once  the  coefficients  of 
the  combinations  constraints  have  been  computed  they 
should  remain  fixed  throughout  the  remainder  of  that 
cutting  phase.  If  the  initial  box  is  reasonably  small, 
propagation  relative  to  a  fixed  set  of  combined  con¬ 
straints  will  converge  on  an  exact  solution  and  only  a 
single  matrix  inversion  is  needed.  When  the  heuristics 
for  selecting  cuts  fail  to  find  any  viable  cuts  the  proce¬ 
dure  branches. 

The  precise  details  of  Newton’s  procedure  for  selecting 
cuts  can  be  found  in  [4].  Although  these  details  may  be 
required  to  exactly  replicate  Newton’s  behavior,  it  seems 
likely  that  the  qualitative  performance  of  the  system  can 
be  duplicated  using  only  the  heuristics  outlined  here. 

4  Conclusion 

This  paper  describes  a  branch  and  propagate  algorithm 
to  find  all  isolated  solutions  for  a  system  of  polynomial 
equations  over  the  reals.  Our  techniques  were  originally 
inspired  by  local  constraint  propagation  techniques  used 
in  AI  and  logic  programming.  The  mathematics  behind 
our  techniques  are  quite  straightforward  and  far  less  so 
phisticated  than  that  underlying  continuation  methods. 
Yet  our  techniques  seem  to  compare  well  with  continua¬ 
tion  methods  on  a  wide  variety  of  benchmarks.  It  seems 
possible  that  branch  and  propagate  techniques  provide 
the  most  effective  approach  to  nonlinear  root  finding  just 
as  they  provide  the  most  effective  approach  to  a  wide  va¬ 
riety  of  other  NP  complete  problems. 
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